!-------------------------------------------------------------------------------------
subroutine set_filter(filter)
!-------------------------------------------------------------------------------------

  use dimensions

  use parameters, only :  kFilterAngle

  implicit none

  real,    intent(out) :: filter(jz,jy,jx)
  real                 :: kk, tantheta
  integer              :: ikx, iky, ikz

      filter=1.0
      tantheta = tan(kFilterAngle)

      do ikx=1,jx
         do iky=1,jy
            do ikz=1,jz
               kk = (kXvec(ikx))**2 + (kYvec(iky))**2 
               !if( kk > kFilter**2 ) filter(ikz,iky,ikx) = 0   !cylindrical
               if( kk > (kZvec(ikz)*tantheta)**2 ) filter(ikz,iky,ikx) = 0   !conical
            enddo
         end do
      end do

end subroutine set_filter

!-------------------------------------------------------------------------------------
subroutine filter_F_erase_nonlin(dfu, fInit)
!-------------------------------------------------------------------------------------

  use fft_3d

  use parameters, only:  modeSave

  implicit none

  complex, intent(inout)   :: dfu(jv, jz, jy, jx, jDir)
  complex, intent(in)      :: fInit(jv,jz,jDir)


  complex    :: dummy3D(jz,jy,jx)
  complex    :: fp(jz)
  complex    :: fm(jz)

  integer    :: iv, iDir, px, py, mx, my


  !-- index of the main mode

  px = modeSave + 1
  mx = jx - modeSave + 1
  py = 1
  my = 1


  do iDir = 1,jDir
     do iv=1,jv

        dummy3D(:,:,:) = dfu(iv,:,:,:,iDir)

        call fft_forward_x(dummy3D)
        call fft_forward_y(dummy3D)

        !-- preserve the main mode

        fp(:)=dummy3D(:,py,px)
        fm(:)=dummy3D(:,my,mx)

        dummy3D(:,:,:) = 0;

        !-- restore the zero's mode and x,y-independent equilibrium
 

        dummy3D(:,py,px) = fp(:)
        dummy3D(:,my,mx) = fm(:)
        dummy3D(:,1,1) = fInit(iv,:,iDir)


        call fft_backward_y(dummy3D)
        call fft_backward_x(dummy3D)

        dfu(iv,:,:,:,iDir) = dummy3D(:,:,:)
 
     end do
  end do

end subroutine filter_F_erase_nonlin

!-------------------------------------------------------------------------------------
subroutine filter_U_erase_nonlin(uxMat, uyMat)
!-------------------------------------------------------------------------------------

  use fft_3d
  use parameters, only:  modeSave, ux0, uy0

  implicit none

  complex, intent(inout)   :: uxMat(jz, jy, jx, jDir)
  complex, intent(inout)   :: uyMat(jz, jy, jx, jDir)

  complex    :: dummy3D(jz,jy,jx)
  complex    :: fp(jz)
  complex    :: fm(jz)

  integer    :: iDir, px, py, mx, my


  !-- index of the main mode

  px = modeSave + 1
  mx = jx - modeSave + 1
  py = 1
  my = 1


  !-- x-velocity

  do iDir = 1,jDir
 
        dummy3D(:,:,:) = uxMat(:,:,:,iDir) - ux0(iDir) 

        call fft_forward_x(dummy3D)
        call fft_forward_y(dummy3D)

        !-- preserve the main mode

        fp(:)=dummy3D(:,py,px)
        fm(:)=dummy3D(:,my,mx)

        dummy3D(:,:,:) = 0;

        !-- restore the main mode
 
        dummy3D(:,py,px) = fp(:)
        dummy3D(:,my,mx) = fm(:)


        call fft_backward_y(dummy3D)
        call fft_backward_x(dummy3D)

        uxMat(:,:,:,iDir) = dummy3D(:,:,:) + ux0(iDir) 

  end do


  !-- y-velocity

  do iDir = 1,jDir
 
        dummy3D(:,:,:) = uyMat(:,:,:,iDir) - uy0(iDir) 

        call fft_forward_x(dummy3D)
        call fft_forward_y(dummy3D)

        !-- preserve the main mode

        fp(:)=dummy3D(:,py,px)
        fm(:)=dummy3D(:,my,mx)

        dummy3D(:,:,:) = 0;

        !-- restore the main mode
 
        dummy3D(:,py,px) = fp(:)
        dummy3D(:,my,mx) = fm(:)


        call fft_backward_y(dummy3D)
        call fft_backward_x(dummy3D)

        uyMat(:,:,:,iDir) = dummy3D(:,:,:) + uy0(iDir) 

  end do

end subroutine filter_U_erase_nonlin

!-------------------------------------------------------------------------------------
 subroutine filter_U_modulation (uxMat, uyMat, filter)
!-------------------------------------------------------------------------------------

  use parameters, only : iuModulation, iunit, Lx, Ly, ux0, uy0
  use fft_3d

  implicit none


  complex, intent(inout)   :: uxMat(jz, jy, jx, jDir)
  complex, intent(inout)   :: uyMat(jz, jy, jx, jDir)

  real,    intent(in)      :: filter(jz,jy,jx)

  complex                  :: du(jz,jy,jx)
  complex, save            :: uxMod(jx)
  complex, save            :: uyMod(jy)          
  logical, save            :: first_call = .true.

  integer                  :: ix, iy, iz, iv, iDir, ikz


  !-- setting modulation vectors --

  if (first_call) then

     if (iuModulation) then
        
        uxMod(:) = exp( - 14 * xvec(:)**6 / (Lx/2)**6 )
        uyMod(:) = exp( - 14 * yvec(:)**6 / (Ly/2)**6 )

        call dealias_1D(uxMod,jx)
        call dealias_1D(uyMod,jy)
 
        if (Lx .eq. 0) uxMod(:) = 1.
        if (Ly .eq. 0) uyMod(:) = 1.
 
     else

        uxMod(:) = 1.
        uyMod(:) = 1.

     endif

     first_call = .false.

  end if


  !-------------------------------------------


  do iDir = 1,jDir

     !-- filtering x-velocity

     du = uxMat(:,:,:,iDir) - ux0(iDir) 


     do ix=1,jx
        du(:,:,ix) = du(:,:,ix)*uxMod(ix)
     end do

     do iy=1,jy
        du(:,iy,:) = du(:,iy,:)*uyMod(iy)
     end do


     call fft_forward_x(du)
     call fft_forward_y(du)

     du(:,:,:) = du(:,:,:)*filter(:,:,:)

     call fft_backward_x(du)
     call fft_backward_y(du)

     uxMat(:,:,:,iDir) = du +  ux0(iDir)


     !-- filtering y-velocity

     du = uyMat(:,:,:,iDir) - uy0(iDir) 

     do ix=1,jx
        du(:,:,ix) = du(:,:,ix)*uxMod(ix)
     end do

     do iy=1,jy
        du(:,iy,:) = du(:,iy,:)*uyMod(iy)
     end do

     call fft_forward_x(du)
     call fft_forward_y(du)

     du(:,:,:) = du(:,:,:)*filter(:,:,:)

     call fft_backward_x(du)
     call fft_backward_y(du)

     uyMat(:,:,:,iDir) = du +  uy0(iDir)

  enddo


  if (jx .eq. 1) uxMat = 0.
  if (jy .eq. 1) uyMat = 0.

end subroutine filter_U_modulation
 

!-------------------------------------------------------------------------------------
subroutine filter_U_lin (uxMat, filter)
!-------------------------------------------------------------------------------------

  use dimensions
  use fft_wrap

  implicit none

  complex, intent(inout) :: uxMat(jz, jy, jx, jDir)
  real,    intent(in)   :: filter(jz,jy,jx)
  complex               :: dummy3D(jz, jy, jx)

  integer               :: iDir


  do iDir=1,jDir

     dummy3D=uxMat(:,:,:,iDir)

     call dfftw_execute_dft(plan_xyzforward, dummy3D, dummy3D)
     dummy3D = dummy3D/(jx*jy*jz)

     dummy3D(:,:,:) = Dummy3D(:,:,:)*filter(:,:,:)
     dummy3D(1,:,:) = 0.0

     call dfftw_execute_dft(plan_xyzbackward, dummy3D, dummy3D)

     uxMat(:,:,:,iDir) = dummy3D

  end do
  
end subroutine filter_U_lin

!-------------------------------------------------------------------------------------
subroutine dealias_1D(f,j)
!-------------------------------------------------------------------------------------

  use parameters, only : iDealias
  use fft_wrap

  implicit none

  integer j
  complex f(j)

  if (.not. iDealias) return

  if (j < 6) return

  if (j .eq. jx) then
     call dfftw_execute_dft(plan_xforward, f, f)
     f=f/float(j)
     f(j/2-j/6 : j/2+2+j/6) = (0.0,0.0)
     call dfftw_execute_dft(plan_xbackward, f, f)
     return
  end if

  if (j .eq. jy) then
     call dfftw_execute_dft(plan_yforward, f, f)
     f=f/float(j)
     f(j/2-j/6 : j/2+2+j/6) = (0.0,0.0)
     call dfftw_execute_dft(plan_ybackward, f, f)
     return
  end if

  if (j .eq. jz) then
     call dfftw_execute_dft(plan_zforward, f, f)
     f=f/float(j)
     f(j/2-j/6 : j/2+2+j/6) = (0.0,0.0)
     call dfftw_execute_dft(plan_zbackward, f, f)
     return
  end if

  return

end


!-------------------------------------------------------------------------------------
subroutine dealias_5D(dfu)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only : iDealias
  use fft_wrap

  implicit none

  complex, intent(inout)   :: dfu(jv, jz, jy, jx, jDir)

  if (.not. iDealias) return 

  if (jx > 5) then 
     call dfftw_execute_dft(plan_x5forward, dfu, dfu)
     dfu = dfu/jx
     dfu( :, :, :, jx/2-jx/6 : jx/2+2+jx/6, :) = 0 
     call dfftw_execute_dft(plan_x5backward, dfu, dfu)
  end if

  if (jy > 5) then
     call dfftw_execute_dft(plan_y5forward, dfu, dfu)
     dfu = dfu/jy
     dfu( :, :, jy/2-jy/6 : jy/2+2+jy/6, :, :) = 0 
     call dfftw_execute_dft(plan_y5backward, dfu, dfu)
  end if

  if (jz > 5) then
     call dfftw_execute_dft(plan_z5forward, dfu, dfu)
     dfu = dfu/jz
     dfu( :, jz/2-jz/6 : jz/2+2+jz/6, :, :, :) = 0 
     call dfftw_execute_dft(plan_z5backward, dfu, dfu)
  end if

  if (jv > 5) then 
     call dfftw_execute_dft(plan_v5forward, dfu, dfu)
     dfu = dfu/jv
     dfu( jv/2-jv/6 : jv/2+2+jv/6, :, :, :, :) = 0
     call dfftw_execute_dft(plan_v5backward, dfu, dfu)
  end if

end

!-------------------------------------------------------------------------------------
subroutine dealias_4D(u)
!-------------------------------------------------------------------------------------

  use dimensions
  use parameters, only : iDealias
  use fft_wrap

  implicit none

  complex, intent(inout)   :: u(jz, jy, jx, jDir)

  if (.not. iDealias) return 

  if (jx > 5) then 
     call dfftw_execute_dft(plan_x4forward, u, u)
     u = u/jx
     u( :, :, jx/2-jx/6 : jx/2+2+jx/6, :) = 0 
     call dfftw_execute_dft(plan_x4backward, u, u)
  end if

  if (jy > 5) then
     call dfftw_execute_dft(plan_y4forward, u, u)
     u = u/jy
     u( :, jy/2-jy/6 : jy/2+2+jy/6, :, :) = 0 
     call dfftw_execute_dft(plan_y4backward, u, u)
  end if

  if (jz > 5) then
     call dfftw_execute_dft(plan_z4forward, u, u)
     u = u/jz
     u( jz/2-jz/6 : jz/2+2+jz/6, :, :, :) = 0 
     call dfftw_execute_dft(plan_z4backward, u, u)
  end if

end


!-------------------------------------------------------------------------------------
